arXiv:1504.07599v3 [math.NA] 23 Mar 2016 


Explicit Strong Stability Preserving Mnltistage Two-Derivative 

Time-Stepping Schemes 


Andrew J. Christlieb^, Sigal Gottlieb^, Zachary Grant^* David G. Seah 

^Department of Computational Mathematics Science and Engineering, Department of Electrical Engineering, 
and Department of Mathematics, Michigan State University 
^Department of Mathematics, University of Massachusetts, Dartmouth 
^Department of Mathematics, U.S. Naval Academy. 

Abstract 

High order strong stability preserving (SSP) time discretizations are advantageous for use with spatial dis¬ 
cretizations with nonlinear stability properties for the solution of hyperbolic PDEs. The search for high order 
strong stability time-stepping methods with large allowable strong stability time-step has been an active area of 
research over the last two decades. Recently, multiderivative time-stepping methods have been implemented with 
hyperbolic PDEs. In this work we describe sufficient conditions for a two-derivative multistage method to be SSP, 
and find some optimal SSP multistage two-derivative methods. While explicit SSP Runge-Kutta methods exist 
only up to fourth order, we show that this order barrier is broken for explicit multi-stage two-derivative methods 
by designing a three stage fifth order SSP method. These methods are tested on simple scalar PDEs to verify the 
order of convergence, and demonstrate the need for the SSP condition and the sharpness of the SSP time-step in 
many cases. 

1 Introduction 

1.1 SSP methods 

When numerically approximating the solution to a hyperbolic conservation law of the form 

G + /(t/). =0, (1) 

difficulties arise when the exact solution develops sharp gradients or discontinuities. Significant effort has been 
expended on developing spatial discretizations that can handle discontinuities [7], especially for high-order methods. 
These discretizations have special nonlinear non-inner-product stability properties, such as total variation stability 
or positivity, which ensure that when the semi-discretized equation 

Ut = F{u), (2) 

(where w is a vector of approximations to U) is evolved using a forward Euler method, the numerical solution 
satisfies the desired strong stability property, 

||u" -1- AtF{u^)\\ < ||u"||, 0 < At < AtpE, (3) 

where || • || is any desired norm, semi-norm, or convex functional. 

In place of the first order time discretization (3), we typically require a higher-order time integrator, but we still 
wish to ensure that the strong stability property < ||m"|| is satisfied, perhaps under a modified time-step 

restriction, where vF is a discrete approximation to U at time t". In [32] it was observed that some Runge- 
Kutta methods can be decomposed into convex combinations of forward Euler steps, so that any convex functional 
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property satisfied by (3) will be preserved by these higher-order time discretizations. For example, the s-stage 
explicit Runge-Kutta method [33], 
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can be rewritten as convex combination of forward Euler steps of the form (3). If all the coefficients aij and /3ij 
are non-negative, and provided aij is zero only if its corresponding is zero, then each stage is bounded by 
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Noting that each -f for < Atps, and by consistency X]j=o ®i,i = li '''® have 

||w"hi|| < ||ri"|| as long as 


At < CAtpB Vi,j, (5) 

where C = min (We employ the convention that if any of the P's are equal to zero, the corresponding ratios 

^■1,3 

are considered infinite.) The resulting time-step restriction is a combination of two distinct factors: (1) the term 
AtpB that depends on the spatial discretization, and (2) the SSP coefficient C that depends only on the time- 
discretization. Any method that admits such a decomposition with C > 0 is called a strong stability preserving 
(SSP) method. 

This convex combination decomposition was used in the development of second and third order explicit Runge- 
Kutta methods [33] and later of fourth order methods [34, 16[ that guarantee the strong stability properties of 
any spatial discretization, provided only that these properties are satisfied when using the forward Euler (first 
derivative) condition in (3). Additionally, the convex combination approach also guarantees that the intermediate 
stages in a Runge-Kutta method satisfy the strong stability property as well. 

The convex combination approach clearly provides a sufficient condition for preservation of strong stability. 
Moreover, it has also be shown that this condition is necessary [4, 5, 11, 12[. Much research on SSP methods 
focuses on finding high-order time discretizations with the largest allowable time-step At < CAtpE by maximizing 
the SSP coefficient C of the method. It has been shown that explicit Runge-Kutta methods with positive SSP 
coefficient cannot be more than fourth-order accurate [20, 28[; this led to the study of other classes of explicit SSP 
methods, such as methods with multiple steps. Explicit multistep SSP methods of order p > 4 do exist, but have 
severely restricted time-step requirements [7[. Explicit multistep multistage methods that are SSP and have order 
p > 4 have been developed as well [17, 2[. 

Recently, multi-stage multiderivative methods have been proposed for use with hyperbolic PDEs [29, 37]. 
The question then arises as to whether these methods can be strong stability preserving as well. Nguyen-Ba 
and colleagues studied the SSP properties of the Hermite-Birkoff-Taylor methods with a set of simplified base 
conditions in [23]. In this work we consider multistage two-derivative methods and develop sufficient conditions 
for strong stability preservation for these methods, and we show that explicit SSP methods within this class can 
break this well-known order barrier for explicit Runge-Kutta methods. Numerical results demonstrate that the 
SSP condition is useful in preserving the nonlinear stability properties of the underlying spatial discretization and 
that the allowable time-step predicted by the SSP theory we developed is sharp in many cases. 

1.2 Multistage multiderivative methods 

To increase the possible order of any method, we can use more steps (e.g. linear multistep methods), more stages 
(e.g. Runge-Kutta methods), or more derivatives (Taylor series methods). It is also possible to combine these 
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approaches to obtain methods with multiple steps, stages, and derivatives. Multistage multiderivative integration 
methods were first considered in [24, 38, 35], and multiderivative time integrators for ordinary differential equations 
have been developed in [30, 31, 14, 15, 22, 25, 3[, but only recently have these methods been explored for use with 
partial differential equations (PDEs) [29, 37[. In this work, we consider explicit multistage two-derivative time 
integrators as applied to the numerical solution of hyperbolic conservation laws. 

We consider the system of ODEs (2) resulting from the spatial discretization of a hyperbolic PDE of the form 
(1). We define the one-stage, two-derivative building block method = u” + aAtF{u") + pAt^F{u") where 
a > 0 and /3 > 0 are coefficients chosen to ensure the desired order. This method can be at most second order, 
with coefficients a = 1 and /3 = |. This is the second-order Taylor series method. To obtain higher order explicit 
methods, we can add more stages: 

i—1 

u’^ + Aty^^ (aijF{y^^^) + AtajjFjy^^^)) , i = l,...,s (6) 

s 

+ Atyy^ {hjF{y^^^) -f AtbjF{y^^^)j . 

We can write the coefficients in matrix vector form, where 

\ 


/ 

We let c = Ae and c = Ae, where e is a vector of ones. These coefficients are then selected to attain the desired 
order, based on the order conditions written in Table 1 as described in [3, 6[. 

Remark 1. In this work, we focus on multistage two-derivative methods as time integrators for use with hyperbolic 
PDEs. In this setting, the operator F is obtained by a spatial discretization of the term Ut = —f{U)x to obtain the 
system ut = F{u). This is the typical method-of-lines approach, and SSP methods were introduced in the context 
of this approach. The computation of the second derivative term F should follow directly from the definition of F, 
where we compute F = F(u)t = FuUt = FuF. In practice, the calculation of Fu may be computationally prohibitive, 
as for example in the popular WENO method where F has a highly nonlinear dependence on u. 

Instead, we adopt a Lax- Wendroff type approach, where we use the fact that the system of ODEs arises from 
the PDE (1) to replace the time derivatives by the spatial derivatives, and discretize these in space. This approach 
begins with the observation that F{u) = ut = Ut -\- 0{Ax'^) (for some integer m). The term F(u) is typically 
computed using a conservative spatial discretization D^ applied to the flux: 

F{u)^D4-fiu)). 
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Next we approximate 

F{u)t = Utt ^ Uu = -f{UU = {-f{U)t), = {-f'{U)Ut)^ « {-f{u)ut) , 

where a (potentially different) spatial approximation Dj, is used. This means that 

F{u)t = Uu = Uu + 0{Ax") ^F + 0{Ax^) 

(for some integers n and r). 
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Since Ft = F -\- 0{l\x^), we will not necessarily obtain the time order AF when we satisfy the order conditions 
in Table 1, as our temporal order is polluted by spatial errors as well. The concern about the Lax-Wendroff approach 
is that the order conditions in Table 1 are based on the assumption that F = Ft, which is not exactly correct in our 
case, thus introducing additional errors. However, these errors are of order r in space, so that in practice, as long 
as the spatial errors are smaller than the temporal errors, we expect to see the correct order of accuracy in time. 

To verify this, in Section 4.2 we perform numerical convergence studies of these temporal methods where F is 
approximated by high order spatial schemes and compare the errors from these methods to those from well-known 
SSP Runge-Kutta methods. We observe that if the spatial and temporal grids are refined together, the expected 
order of accuracy is demonstrated. Furthermore, if the spatial grid is held fixed but the spatial discretization is 
highly accurate, the correct time-order is observed until the time error falls below the spatial error. We conclude 
that in practice, it is not necessary to compute Ft exactly, and that the use of a Lax-Wendroff type procedure that 
replaces the temporal derivatives by spatial derivatives and discretizes each of these independently, does not destroy 
the temporal accuracy. 

This observation is not new: many other methods have adopted this type of approach and obtained genuine 
high-order accuracy. Such methods include the original Lax-Wendroff method [21], ENO methods [9], finite volume 
ADER methods [36], the finite difference WENO Schemes in [27], and the Lax-Wendroff discontinuous Galerkin 
schemes ]26]. 

In the next section we will discuss how to ensure that a multistage two-derivative method will preserve these 
strong stability properties. 


p = 1 

= 1 

p = 2 

-1- 6^e = 1 

CO 

II 

b^c^ -h 2b^c = 1 

-1- b"'"c -\- c = 1 

II 

b^c^ A- 3b^c^ = 1 

b"'^cAc -1- 6^cc -1- b^c^ -\- iFAc -\- c = | 

-h 26^ic -t b^c^ = ^ 
b'^A^c -\- Ac -1- -1- P" Ac = T 

II 

-h 4b'^c^ = 1 

b^c^Ac A- b^c^c A- -|- 2b^cAc + 2b^cc = T 

b^cAc^ -A 2b^ cAc + lAc^ -1- lA Ac^ -A 2bA Ac = ^ 

b^cA^c -A b^cAc + b^cAc -1- iFcAc -A iFcc -A iF A^c -1- bF Ac -1- hFAc = ^ 

b'^(Ac)(Ac) -A 2b^cAc -A F + 2lFcAc + 2b^cc = ^ 

b'^Ac^ -A 3b^Ac^ -A iFc^ = ^ 

b^A(cAc) -A b^ A(cc) + b^Ac^ -A b^AAc -A b^Ac -A iFcAc -A iFcc = T 
b'^A^c^ -1- 2b'^ AAc + b^Ac^ -A iF Ac^ -1- 2iFAc = T 

-1- b^ A^c + b^ AAc + b^ AAc -1- b'^ Ac + iF A^c -A iF Ac -A iF Ac = ^ 


Table 1: Order conditions for multi-stage multiderivative methods of the form (6) as in [3]. 
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2 The SSP condition for multiderivative methods 


2.1 Motivating Examples 

To understand the strong stability condition for multiderivative methods, we consider the strong stability properties 
of a multiderivative building block of the form 

= y- + aAtF{u”) + PAt^Fiu^), 

and begin with the simple linear one-way wave equation Ut = U^- This equation has the property that its second 
derivative in time is, with the assumption of sufficient smoothness, also the second derivative in space: 


Utt = {U^)t = {Ut). = u. 


We will use this convenient fact in a Lax-Wendroff type approach to define F{vF) by a spatial discretization of 

UXX • 

For this problem, we define F by the original first-order upwind method 


F{u")j ■- ^ «+i - m") « Ux{xj), 
and F by the second order centered discretization to U..- 

F{u")j ■- ^ «+i - 2u7 -b u"_i) « Uxx{xj). 

These spatial discretizations are total variation diminishing (TVD) in the following sense: 

= u" + AtF{u") is TVD for At < Ax, 


n + 1 


= u" -b Ai?F{u^) is TVD for At < 


(7a) 

(7b) 

(8a) 

(8b) 


Remark 2. Note that we chose a second derivative F in space that is not an exact derivative Ft{u") of F{u'^) in 
the method-of-lines formulation. However, as noted in Remark 1 and will be shown in the convergence studies, if the 
spatial and temporal grids are co-refined, this provides a sufficiently accurate approximation to F{u^). The exact 
derivative can be obtained by applying the upwind differentiation operator to the solution twice, which produces 

Ft ■- (Mi+2 - 2u"+i -b m") . 

However, computing F = Ft using this formulation does not satisfy the condition (8b) for any value of At. 

To establish the TVD properties of the multiderivative building block we decompose it: 

= u" + aAtF{u^) -b pAt^F{u") — au" -b aAtF{u^) -b (1 — a)u'^ + PAt^F{u^) 

= a(u" + ^AtF{u’^)'j -b (1 - a) -b j^At^F{u^)'j . 

It follows that for any 0 < o < 1 this is a convex combination of terms of the form (8a) and (8b), and so 

«AfF’(u"))|| -b(l-a) fu"-b T^Af"F’(u")^ 

IV a / Wtv \ 1 — a J 


II 71 + 1 II 

||ii \\tv S ^ 


_ II nil I /i \ II FLU _ II nil 

< a ||u -b (1 — a) ||u lljn^ < ||m ||tv 
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for time-steps satisfying At < ^Aa; and At^ < Aa;^. The first restriction relaxes as a increases while the second 

becomes tighter as a increases, so that the value of a that maximizes these conditions occurs when these are equal. 
This is given by 


2 2 
2 a 


a^ycfi -|- 8/3 — c? 
4^ 


Using this SSP analysis, we conclude that 


-I- alS.tF{u^) /3At^T’(u") 


^ II "II f ^ 

< \\u for At < — ---r- Ax. 

- II IITl^ - 4^ 


(9) 


Of course, for this simple example we can directly compute the value of At for which the multiderivative building 

At 


block is TVD. That is, with Sa- — we observe that 


II n+1 II 

\\u \\tv 


provided that 


(1 — aA — 2/3A^)Mj -I- (aA -|- /3A^)u"_,_i + /3A^uJ'_i) 

y/a'^ + SP — a 


Ity < l|w \\tv, 


l-aX- 2py > 0 


A < 


4/3 


We see that for this case, the SSP bound is sharp: the convex combination approach provides us exactly the same 
bound as directly computing the requirements for total variation. 

We wish to generalize this for cases in which the second derivative condition (8b) holds for At < KAIfe where 
K can take on any positive value, not just For the two-derivative building block method this can be done 
quite easily: 

Theorem 1. Given F and F such that 


and 


the two-derivative building block 


\\u" + AtF{u")\\ < ||«” 
\\u” + AyF{u")\\ < ||u" 


for At < AtpB, 


for At < KAtpE, 


n-l-l 


= u" -b aAtF{u") + PAyF{u^) 


satisfies the monotonicity condition ||u"+^|| < Hu"'!! under the time-step restriction 

(ya'^K'^ + 4/3 - AtpB. 




Proof. As above, we rewrite 


, y -b ytF{uy -b (1 - a) -b r^AyPiu’-^)^ , 


which is a convex combination provided that 0 < a < 1. The time-step restriction that follows from this convex 
combination must satisfy ^At < Atps and j^AP < K^AtpF- The first condition becomes At < f^AtFF while 

the second is At < ^^-f^KAtFE- We observe that on 0 < a < 1 the first term encourages a larger a while the 
second term is less restrictive with a smaller a. The two conditions are balanced, and thus the allowable time-step 
is maximized, when we equate the right hand sides: 


a 

a 


1 — a 


K 


a = ^ ya^K^ + Ap - aX) . 


Now using the first condition. At < -Atpp we obtain our result. 


□ 
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A more realistic motivating example is the unique two-stage fourth order method 


At Af^ ■ 

* n I T-t/ n\ , TTi/ n\ 

u = u + )> 

+AtJ^(u") +^(J’(u") + 2F(u’)). (10) 

The first stage of the method is a Taylor series method with ^, while the second stage can be written 

/\+2 

„n+i ^ +AtT'(n")-b ^(F(M")-b2F(M*)) 

= a + (1 - + 2 F(m‘)) 

= (1 - a) + ^^AtF(n") + At"F(n")^ + a + ^F{u*)j . 


For 0 < a < 1 this is a convex combination of two terms. The first term is of the form (9), which gives the time-step 
restriction 


At < 


6 


4- 3a 



AtpE- 


The second is of the form (8b), so we have At < y ^ AtpE- We plot these two in Figure 1, and we observe that the 
first term is decreasing in a (blue line) while the second term is increasing in a (red line). As a result, we obtain the 
optimal allowable time-step by setting these two equal, which yields a « 0.3072182638002141 and the corresponding 
SSP coefficient, C ~ 0.6788426884782078. A direct computation of the TVD time-step for this case, which takes 
advantage of the linearity of the problem and the spatial discretization, gives the bound At < (-^/S — 1) > 0.6788. 
This shows that, as we expect, the SSP condition is not always sharp. 

The convex combination approach becomes more complicated when dealing with multi-stage methods. It is the 
most appropriate approach for developing an understanding of the strong stability property of a given method. 
However, it is not computationally efficient for finding optimal SSP methods. In the following section, we show 
how to generalize the convex combination decomposition approach and we use this generalization to formulate an 
SSP optimization problem along the lines of [16, 19]. 


2.2 Formulating the SSP optimization problem 

As above, we begin with the hyperbolic conservation law (1) and adopt a spatial discretization so that we have the 
system of ODEs (2). The spatial discretization F is specially designed so that it satisfies the forward Euler (first 
derivative) condition 

Forward Euler condition ||u" + AtF{u^)\\ < ||m"|| for At < AtpE, (11) 

for the desired stability property indicated by the convex functional || ■ ||. For multiderivative methods, in addition 
to the first derivative, we need to appropriately approximate the second derivative in time Utt, to which we represent 
the discretization as F. It is not immediately obvious what should be the form of a condition that would account 
for the effect of the At^F term. Motivated by the examples in the previous sections, we choose the 

Second derivative condition ||u”'-|-At^F’(M")|| < ||it”'|| for At < KAtpE, (12) 

where K is a scaling factor that compares the stability condition of the second derivative term to that of the forward 
Euler term. Given conditions (11) and (12), we wish to formulate sufficient conditions so that the multiderivative 
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method (6) satisfies the desired monotonicity condition under a given time-step. First, we write the method (6) in 
an equivalent matrix-vector form 

y = en" + AtSF(y) + At^SFiy), (13) 

where 


'A 0 


A 0 

. 0 . 

and S = 

o-> ' 
H ' 

O 1 

1 _ 


and e is a vector of ones. We are now ready to state our result: 

Theorem 2. Given spatial discretizations F and F that satisfy (11) and (12), a two-derivative multistage method 
of the form (13) preserves the strong stability property ||u"^^|| < ||rt’^|| under the time-step restriction At < rAtpE 
if satisfies the conditions 

e > 0 (14a) 

5 > 0 (14b) 

S>0 (14c) 

for some r > 0. In the above conditions, the inegualities are understood component-wise. 

Proof. We begin with the method (13), and add the terms rSy and rSy to both sides to obtain 

(^I-\-rS-IrS'^y = u”e + rS (y-\-^F{y)^ -\-fS ^^F{y)^ , 
y = i?(eu") -I- P ^y -I- ^F’(y)^ -I- g ^y -I- ^F(y)^ , 

where 

R= (^I-\-rS-\-fSy\ P^rRS, Q = fRS. 

If the elements of P, Q, and Re are all non-negative, and if i? -|- P -|- Q = 7, then these three terms describe a 
convex combination of terms which are SSP, and the resulting value is SSP as well 

At AF . 

liyil < R\\eu-\\ + P\\y + —P(y)|i + g|ly + ^P(y)||, 

under the time-step restrictions At < rAtpE and At < Ky/fAtYE. As we observed above, the optimal time-step 
is given when these two are set equal, so we require r — Ky/f. Conditions (14a)-(14c) now ensure that P > 0, 
Q > 0, and Re > 0 component-wise for f = and our method (13) preserves the strong stability condition 
||^n+i|| ^ ||u"|| under the time-step restriction At < rAtps. □ 

This theorem gives us the conditions for the method (13) to be SSP for any the time-step At < rAtpE. This 
allows us to formulate the search for optimal SSP two-derivative methods as an optimization problem, similar 
to [16, 18, 7], where the aim is to find C — maxr such that the relevant order conditions (from Section 1.2) 
and SSP conditions (14a)-(14b) are all satisfied. Based on this, we wrote a Matlab optimization code for finding 
optimal two-derivative multistage methods [8], formulated along the lines of David Ketcheson’s code [19] for finding 
optimal SSP multistage multistep methods in [17, 2[. We used this to find optimal SSP multistage two-derivative 
methods of order up to p = 5. However, we also used our observations on the resulting methods to formulate closed 
form representations of the optimal SSP multistage two-derivative methods. We present both the numerical and 
closed-form optimal methods in the following section. 








Remark 3. An alternative approach to defining a MSMD SSP method is to begin with spatial discretization that 
satisfies the “Taylor series” condition, defined by 

\\u”' + AtF{u”') + ^At'^F{u”)\\ < ||ii"|| for At < AtTS = K 2 AtFE- (15) 

This condition replaces (12), and allows us to rewrite (6) as convex combinations of forward Euler and Taylor 
series steps of the form (15). A similar optimization problem can he defined based on this condition. This approach 
was adopted in [23]. 

This condition is more restrictive than what we consider in the present work. Indeed, if a spatial discretization 
satisfies conditions (11) as well as (12), then it will also satisfy condition (15), with K 2 = + 2 — K. However, 

some methods of interest cannot be written using a Taylor series decomposition, including the two-stage fourth 
order method in (10). For these reasons, we do not explore this approach further, but we point out that we have 
experimented with this alternative formulation to generate some optimal SSP methods. Henceforth, we restrict 
our attention to spatial discretizations that satisfy (11) and (12). 


3 Optimal SSP multiderivative methods 

3.1 Second order methods 

Although we do not wish to use second order methods in our computations, it is interesting to consider the strong 
stability properties of these methods both as building blocks of higher order methods and as simple example that 
admit optimal formulations with simple formulas. 

One stage methods. The second-order Taylor series method, 

u”'+^ =u”’+ AtF{u”') + ]^At'^F[u”‘), (16) 

is the unique one-stage second-order method. Theorem 1 for the building block (9) with coefficients a = 1 and 
/3 = I gives the condition At < CAtpE with C = K\/K^ + 2 — K^. Unlike the SSP single derivative Runge-Kutta 
methods, the SSP coefficient is not just dependent on the time-stepping method, but also on the value K which 
comes from the second derivative condition (12). As noted above, the Taylor series method can serve as the 
basic building block for two-derivative methods, just as the first-order FE can be used to build higher-order RK 
integrators. 

Two stage methods. The optimal SSP two stage second order methods depends on the value of K in (12). A 
straightforward SSP analysis via convex combinations, and solving the equations for the order conditions allows 
us to formulate the optimal methods without using the optimization code [8]. However, we used the optimization 
code to verify our results. 

\i K < the optimal methods are of the form 

It* = m" -I- -AtF{u”), 
r 

= M" + iAt(F(u")-fF(u*))-b^^At2F(M"), (17) 

z Zr 

where 

r=^(l-K'^ + Vl + 67^2 . 

These methods are SSP for C = r, where r > 1 whenever K > 0. Notice that if A" = 0 we have r = 1 and our 
method reduces to the standard two stage second order Runge-Kutta method. 
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As K increases, the time-step restriction imposed by the condition (12) is alleviated, and consequently the 
optimal method includes more of the second derivative terms. If A > < /1 we get an optimal method 


Z o 

u F \mF{u) F \l^t^F(u), (18) 

2 8 

which is simply two Taylor series steps with | At in each one, and so has C = 2Ky/K^ -I- 2 — 2K^. We see in Figure 
2 the SSP coefficient C vs. K for the two methods above (in green and red respectively), and for many numerically 
optimal methods In blue. 


u = 

n + 1 _ 



a 

Figure 1: The time-step restriction 
for the two-stage fourth order method 
is the minimum of two terms, one of 
which is decreasing while the other is 
increasing in a. 



Figure 2: The SSP coefficient C as a 
function of the coefficient K in (12) for 
the optimal two-stage two-derivative 
second order methods given in Section 
3.1. 


Note that (18) has four function evaluations compared to the three function evaluations above in (17). If we 
assume all the function evaluations cost the same, it will never pay off to use the second method, as the first 
method is always more efficient. However, if one has a special case where the cost of computing F is negligible the 
second method may still be worthwhile. 

For this method and all others, once we have the optimal Butcher arrays and the SSP coefficient C = r, we can 
easily convert them to the Shu-Osher form as follows: 

y. Given A, Ahat, b, bhat, and r, the Shu-Osher matrices are given by 
z=0.0*b; I=eye(3);e=ones(3,1); 

S= [A, z; b ’ , 0] ; Shat= [Ahat, z; bhat ’ , 0] ; 

Ri=(I+r*S+(r-2/k-2)*Shat); 
v=Ri\e; 

P = r*(Ri\S); 

Q= r-2/k-2*(Ri\Shat); 
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3.2 Third order methods 


Many two-stage two-derivative third order methods exist. As above, the optimal C depends on the value of K 
in (12). Through numerical search using our optimization code [8], we found optimal methods for the range 
0.1 < A < 5. Figure 3 (blue line) shows the SSP coefficient C vs. K. We also solved the order conditions explicitly 
and analyzed the SSP conditions (14a)-(14c) to Hnd the values of the coefficients of these methods as functions of 
r = C and K. These methods all have the form 


u* = u’^ + aAtF{u”') + aAt^ F{u"), 

= u" + biAtF{u^) + b 2 AtF{u*) + biAt^ F{u") + b 2 At^ F{u*), 


where the coefficients satisfy 


a = I {KVIO + 2 - K^) 

^ K\/ K^+2+K'^ 3if2 

bl = I - ^ab2 - ^ 


a = 


1 

2 


a 


2 


6l = 1 — &2 


62 


1 

6a 


1 

2 


ab2 


(19) 


( 20 ) 


Note that the first stage is always a Taylor series method with At replaced by aAt, which is underscored by 
the fact that SSP coefficient is given by r = ^ -I- 2 — The SSP conditions (14a)-(14c) provide the 

restrictions on the size of r, and thus on C = r. We observe that the condition that restricts us most here is the 
non-negativity of Qs.i, and so we must select a value of r such that this value is zero. To do this, we select r to be 
the smallest positive root of po -|- pir + p 2 r^ -I- pzr^ where the coefficients pi are all functions of K 


Po 

Pi 


2A'(VA'2 -f 2 - 3K) + 4A'^(VA'2 -h2 - K) 

-ao, P2 = (1 - ao)/(2A'^), Pa 


^ + K 


6K3 


Some examples of the value of r as a function of K are given in Table 2. The Matlab script for this method is 
found in Appendix B. 


K 0.25 

0.4 

0.5 

0.6 

0.7 

0.8 

1.0 

1.25 

1.5 

1.75 

2.5 

3 

3.5 

4 

r 0.48 

0.71 

0.84 

0.94 

1.03 

1.11 

1.23 

1.33 

1.39 

1.44 

1.51 

1.54 

1.55 

1.56 


Table 2: The SSP coefficient C = r for each K in the two-stage third order method (19). 


The Shu-Osher representations of the methods can be easily obtained by using r and K to solve for f, R, P, 
and Q. For example, for K — ^ we have r = 1.04 and the method in Butcher form has 


a = 0.594223212099088, _ 0.693972512991841 

a = 0.176550612898679, ^ [ 0.306027487008159 


0.128597465450411 

0.189553898228989 


The optimal Shu Osher formulation for these values is given by Re = (1,0,0)^, 


P = 

Q = 


0 

0.618033988749895 

0.271611333775367 

0 

0.381966011250105 

0 


0 

0 

0.318290138472780 

0 

0 

0.410098527751853 


0 

0 

0 

0 

0 

0 
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Figure 3: The SSP coefficient vs. K for 
two stage methods. Third order methods 
are in blue, the fourth order methods are 
red. 



Figure 4: The SSP coefficient vs. K for 
three stage methods. Fourth order meth¬ 
ods are in blue (top line), the fifth order 
methods in red (bottom fine). 


3.3 Fourth order methods 


3.3.1 Fourth order methods: The two-stage fourth-order method 


The two-stage two-derivative fourth order method (10) is unique; there is only one set of coefficients that satisfy 
the fourth order conditions for this number of stages and derivatives. The method is 


At At^ ■ 

* n I 771 / n\ , 7-1/ n\ 

u =u +^F{u ) + —F{u ) 

= u" -b AtF{u’^) + ^(F(u") -b 2F{u)). 


The first stage of the method is a Taylor series method with ^, while the second stage can be written as a linear 
combination of a forward Euler and a second-derivative term (but not a Taylor series term). The SSP coefficient 
of this method is larger as K increases, as can be seen in Figure 3. 

To ensure that the SSP conditions (14a)-(14c) are satisfied, we need to select the largest r so that all the terms 
are non-negative. We observed from numerical optimization that in the case of the two-stage fourth order method 
the term {Re)^ gives the most restrictive condition: if we choose r to ensure that this term is non-negative, all the 
other conditions are satisfied. Satisfying this condition, the SSP coefficient C = r is given by the smallest positive 
root of the polynomial: 

(i?e)3 = -b dlfV® - 12A'V^ - 24E"‘r -b 24,K^. 


The Shu-Osher decomposition for the optimal method corresponding to this value of K is 


= 1 - 


-b 


n + 1 


„2 


At. 


-F{u-) ]+—{u- + t^AFFiu") 




( 21 ) 






-F{un + 


r2(4/s:2 _^2) 

247^4 




u" + t^At^F{u-) )+^(n* + t^At^F{u-) 




3.3.2 Fourth order methods: Three stage methods 

If we increase the number of stages to three, we can construct entire families of methods that obtain fourth-order 
accuracy, and are SSP with a larger allowable time-step. For these methods, we were not able to find closed 
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form solutions, but our optimization code [ 8 ] produced methods for various values of K. The SSP coefficient as 
a function of K for these methods is given in Figure 3, and we give the coefficients for selected methods in both 
Butcher array and Shu-Osher form in Appendix A. 


3.4 Fifth order methods 


As mentioned above, it was shown that explicit SSP Runge-Kutta methods cannot have order p > 4 [20, 28]. This 
order barrier is broken by multiderivative methods. If we allow three stages and two-derivative we can obtain a fifth 
order SSP method. The explicit three-stage fifth order method has twelve coefficients and sixteen order conditions 
that need to be satisfied. This is possible if some of the coefficients are set to zero, which allows several of the 
order conditions to be repetitive and satisfied automatically. The methods resulting from our optimization routine 
all had the simplihed form 


u* = + a 2 iAtF{u’^) + a 2 iAt^ F{u^) 

u** = u" -I- a3iAtF{u’^) + a3\At^F{vF) + a32At^F{u*) ( 22 ) 

= u" + AtF{u'") +Af {biF{u'^) + b2F{u) + h3F{u*)y 

The coefficients of the three-stage fifth order method are then given as a one-parameter system, depending only 
on 021 , that are related through 


U 21 


032 

&2 


1 2 

2 « 21 , 

1 

To 


U 31 


3/5 - 021 
1 - 2021 ’ 

( 5 - 021 ) 5 


— O 2 I 


021(1 — 2 o 2 i )3 (1 — 202 i )^ 


031 


2031 — 1 

12 o 2 i (031 — 021 ) ’ 


63 = 


1 — 2021 


12031 (031 — 021 ) ’ 


1 ( 1 - 021 )^ . 

2 (1 - 202 i )2 

61 = - — 62 — 53 - 


To satisfy the SSP conditions (14a)-(14c), we must ensure that {Re )3 is non-negative. Based on the optimization 
code we observed that the extreme case of {Re )3 = 0 gives the optimal methods, and we can obtain 021 as a function 
of K and r through 


021 = -^ 1 + TTjr + T^r - 


^"' + :^"'-^^'- 240r + 240 


Now, we wish to ensure that Qa^i is nonnegative. The SSP coefficient C = r is then chosen as the largest positive 
root of 


Qs.i = 10r^a2i — (lOOif^ -I- 10r^)a2i -I- (1301^^ -I- 3r^)a2i — 50A'^a2i + . 


The Matlab script in Appendix C solves for the largest r that satishes the SSP conditions (14a)-(14c), and then 
computes the coefficients of the optimal methods both in Butcher array and Shu-Osher form. This approach yields 
the same optimal methods as those obtained by our optimization code [ 8 ]. In Figure 5 we show values of 021 and 
r for given values of K. 

Once again, the Shu-Osher decomposition is needed for the method to be SSP, and is easily obtained. For 
example, for K = ^ we have r = 0.6747 and the method becomes 


■ 1.0 


■ 0 

0 

0 

0 ■ 

0.2369970626512336 

, P = 

0.5064804704259125 

0 

0 

0 

0.7810723816004148 


0.1862033791874200 

0 

0 

0 

. 0.0 


. 0.5769733539128722 

0 

0 

0 . 
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K 

021 

C 

K 

021 

C 

0.1 

0.7947 

0.1452 

1.1 

0.7393 

0.8114 

0.2 

0.7842 

0.2722 

1.2 

0.7374 

0.8335 

0.3 

0.7751 

0.3814 

1.3 

0.7359 

0.8523 

0.4 

0.7674 

0.4741 

1.4 

0.7346 

0.8683 

0.5 

0.7609 

0.5520 

1.5 

0.7334 

0.8819 

0.6 

0.7555 

0.6171 

1.6 

0.7324 

0.8937 

0.7 

0.7510 

0.6712 

1.7 

0.7316 

0.9039 

0.8 

0.7472 

0.7162 

1.8 

0.7309 

0.9127 

0.9 

0.7441 

0.7537 

1.9 

0.7302 

0.9205 

1.0 

0.7415 

0.7851 

2.0 

0.7296 

0.9273 


Figure 5: SSP coefficients for three-stage fifth-order methods. Left: the SSP coefficient as a function of K for 
three stage fifth order methods. Right: a table of 021 and the SSP coefficient C, for different values of K as defined 
in (22). The code to generate the coefficients in Butcher and Shu-Osher form is given in the appendix. 


Q = 


0 

0.2565224669228537 

0 

0.0615083849004797 


0 

0 

0.0327242392121651 

0.0803574544380432 


0 0 

0 0 

0 0 

0.2811608067486047 0 


These coefficients as well as the coefficients for the optimal method for any value of K can be easily obtained 
to high precision by the Matlab code in Appendix C. 


4 Numerical Experiments 

4.1 Numerical verification of the SSP properties of these methods 

4.1.1 Example 1: Linear advection with first order TVD spatial discretization 

As a first test case, we consider linear advection, Ut — Ux = 0, with a first order finite difference for the first 
derivative and a second order centered difference for the second derivative defined in (7) 


nu )j := 




Aa; 


Ux{xj), and F{u^)j := 


u'^+i - 2u] + u"_i 


Uxx{Xj ) ■ 


Recall from (8) that this first-order spatial discretization satisfies the 


Forward Euler condition ~ ^ “ ’^1) TVD for At < AIfe = Ax, and the 


Second Derivative condition u 




+ (li)" (“ 1+1 - 2m" + m"-i) is TVD for At < ^AtpE- 


Fo initial conditions, we use a step function 

uo{x) = 


1 if i < 2 ; < |, 
0 otherwise, 


(24) 


with a domain x £ [0,1] and periodic boundary conditions. This simple example is chosen as our experience has 
shown [7] that this problem often demonstrates the sharpness of the SSP time-step. 
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Stages 

Order 

Predicted C 

Observed C 

Stages 

Order 

Predicted C 

Observed C 

1 

2 

0.6180 

0.6180 

2 

4 

0.6788 

0.7320 

2 

2 

1.2807 

1.2807 

3 

4 

1.3927 

1.3927 

2 

3 

1.0400 

1.0400 

3 

5 

0.6746 

0.7136 


Table 3: Comparison of the theoretical and observed SSP coefficients that preserve the nonlinear stability 
properties in Example 1. 


For all of our simulations, we use a fixed grid of size Ax = and a time-step At = XAx where we vary 

0.05 < A < 1. We step each method forward by = 50 time-steps and compare the performance of the various 
time-stepping methods constructed earlier in this work, for K = ^. We test this problem using the two stage third 
order method (19), the two stage fourth order method (10) and the three stage fourth order method in Appendix 
A, the fifth order method (23). We also consider the non-SSP two stage third order method, 

u"-AtF{u") + ^AfF{u"), 

u" - ^AtF{u") + '^AtF{u) + ^At^F{u") + ^At^F{u), (25) 

To measure the effectiveness of these methods, we consider the maximum observed rise in total variation, defined 

by 

max (||u”^^||tv — ||m"||tv) ■ (26) 

We are interested in the time-step in which this rise becomes evident (i.e. well above roundoff error). Another 
measure that we use is the rise in total variation compared to the total variation of the initial solution: 

max (||u”^^||tv — ||w°||tv) ■ (27) 

0<n<]V-l '■ II II II y 

Figure 6 (top) shows the maximal rise in total variation for each CFL value On the left we have the maximal 
per-step rise in TV (26) and on the right, the maximal rise in TV compared to the TV of the initial solution (27). 
We clearly see that once the CFL value passes a certain limit, there is a sharp jump in the total variation of the 
solution. We are interested in the value of ^ at which the time-stepping method no longer maintains the nonlinear 
stability. The fifth order method (black) has the most restrictive value of At before the total variation begins to 
rise. The next most restrictive is the two-stage fourth order method, followed by the two stage third order method. 
The two-stage second order methods have more freedom than these methods, and so have a much larger allowable 
At, while the three stage fourth order, having the most freedom in the choice of coefficients, outperforms all the 
other methods. Figure 6 (bottom) compares the performance of the non-SSP method to the SSP method. This 
graph clearly shows the need for the SSP property of the time-stepping, as the absence of this property results in 
the loss of the TVD property for any time-step. 

We notice that controlling the maximal rise of the total variation compared to the initial condition (27) requires a 
smaller allowable time-step, so we use this condition as our criterion for maximal allowable time-step. A comparison 
of the predicted (i.e. theoretical) values of the SSP coefficient and the observed value for the Taylor series method, 
the two stage methods of order p = 2,3,4, and the three-stage fourth and fifth order methods are shown in Table 3 
We note that for the Taylor series method, the two-stage second order method, the two stage third order method, 
and three stage fourth order method, the observed SSP coefficient matches exactly the theoretical value. On 
the other hand, the two-stage fourth order and the three-stage fifth order, both of which have the smallest SSP 
coefficients (both in theory and practice), have a larger observed SSP coefficient than predicted. For the two-stage 
fourth order case this is expected, as we noted in Section 2.1 that the TVD time-step for this particular case is (as 
we observe here) (^/3 — 1), larger than the more general SSP timestep C — 0.6788. 


u = 

n + 1 _ 
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Figure 6: The rise in total variation as a function of the CFL number. On the left is the maximal per 
time-step rise (26) and on the right the maximal TV rise above the initial TV (27). Top: Comparison of 
a variety of SSP methods. Bottom: Comparison of two-stage third order SSP and non-SSP methods. 


4.1.2 Example 2: MSMD methods with weighted essentially non-oscillatory (WENO) meth¬ 
ods 


The major use of MSMD time-stepping would be in conjunction with high order methods for problems with shocks. 
In this section we consider two scalar problems: the linear advection equation 


Ut -f Ux — 0 


(28) 


and the nonlinear Burgers’ equation 

Ut + = 0 (29) 

on a: G (—1,1). In both cases we use the step function initial conditions (24), and periodic boundaries. We use 
N = 201 points in the domain, so that Ax = 
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Figure 7: Comparison of the rise in total variation as a function of the CFL number for Example 2. Linear 
advection on left and Burgers’ equation on right. The top graphs compare the performance of different 
SSP methods, while the bottom graphs compare the two-stage third order SSP and non-SSP methods. 


We follow our previous work [29] with a minor modification for the spatial discretization. The spatial discretiza¬ 
tion is performed as follows; at each iteration we take the known value u" and compute the flux f{u") = in 
the linear case and f{u"') = | for Burgers’ equation. Now to compute the spatial derivative f{u'^)x we use 

the WENO method [13]. In our test cases, we can avoid flux splitting, as f'{u) is strictly non-negative (below, we 
refer to the WENO method on a flux with f'{u) > 0 as WENO^ and to to the corresponding method on a flux 
with f'{u) < 0 as WENO”). 

Now we have the approximation to Ut at time C, and wish to compute the approximation to Uu- In previous 
work we defined the higher order derivative using central differences, but we have found that additional limiting, in 
the form of the WENO“ differentiation operator, is needed to achieve a pseudo-TVD like property. For the linear 
flux, this is very straightforward as Utt = Uxx- To compute this, we take Ux as computed before, and differentiate 
it using the WENO” method. Now we can compute the building block method. 

For Burgers’ equation, we have Utt = — {UUt)^- We take the approximation to Ut that we obtained above 
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using WENO'^^ we multiply it by vE and differentiate in space using WENO~. The choice of WENO^ followed 
by WENO~ is made by analogy to the first order finite difference for the linear advection case, where we use a 
differentiation operator D"*" followed by the downwind differentiation operator D~ to produce a centered difference 
for the second derivative. The second derivative condition (12) was satisfied by this approach. Now we compute 
the building block method with the approximations to Ut and Utt- In pseudocode, the building block calculation 
takes the form: 

= WENO+ifiu^))-, 

w"; =/'(u")nr u:, = WENO-{f{u")t) 

u" + oAtu" + 

We use the two stage third order SSP method (19), and the non-SSP method (25) the two stage fourth order method 
(10) and the three stage fourth order method in Appendix A, the fifth order method (23). In these simulations, 
we use At = XAx where 0.05 < A < 1.6, and step up to Tjinai = 1.0. At each time-step we compute (27), the 
maximal rise in total variation compared to the total variation of the initial solution. In Figure 7 we observe similar 
behavior to those of the linear advection with first order time-stepping, and once again see that the SSP method 
is needed to preserve the nonlinear stability of WENO as well. 

4.2 Convergence studies 

As a final test case, we Investigate the accuracy of the proposed schemes in conjunction with various high order 
spatial discretization operators. We perform several tests that demonstrate that these methods converge with the 
correct order for linear and nonlinear problems. In the first study (Example 3), we refine the grid only in time, and 
show that if the spatial discretization is sufficiently accurate the multi-derivative methods exhibit the design-order 
of convergence. We also compare the performance of the third order multi-derivative method to the three stage 
third order explicit SSP Runge-Kutta method (SSPRK3,3) [32] and show that the convergence properties are the 
same, indicating that the additional error in approximating Et does not affect the accuracy of the method. In the 
second study (Example 4) we co-refine the spatial and temporal grid by setting At = XAx for a fixed A, and shrink 
Ax. We observed that since the order of the spatial method Is higher than the order of the temporal discretization, 
the time-stepping method achieves the design-order of accuracy both for linear and nonlinear problems. 

Example 3a: temporal grid refinement with pseudospectral approximation of the spatial derivative. 
We begin with a linear advection problem Ut + Ux = 0 with periodic boundary conditions and initial conditions 
uo(x) = 0.5 -I- 0.5sin(a;) on the spatial domain x € [0, 27r]. We discretize the spatial grid with N = 41 equidistant 
points and use the Fourier pseudospectral differentiation matrix D [10] to compute E ~ —Ux ~ —Du. We use 
a Lax-Wendroff approach to approximate Ft ~ Uxx ~ D^u. In this case, the solution is a sine wave, so that 
the pseudospectral method is exact. For this reason, the spatial discretization of F is exact and contributes no 
errors, and the second derivative Ft is also exact Ft = utt = —Dut = —DE = D^u = F. We use a range of time 
steps. At = AAr where we pick A = 0.8, 0.7,0.6,0.5,0.4, 0.3, 0.2, 0.1, and 0.05 to compute the solution to final time 

Tf = 2 . 0 . 

In Table 4 we list the errors for the SSPRK3,3, the two derivative two stage third order method in Section 3.2 
with K = ^ (listed as 2s3p in the table), the unique the two derivative two stage fourth order method (2s4p), 
and the the two derivative three stage fifth order method (22) (3s5p). We observe that the design-order of each 
method is verified. It is interesting to note that the SSPRK3,3 method has larger errors than the multiderivative 
method 2s3p, demonstrating that the additional computation of F does Improve the quality of the solution. 
Example 3b: temporal grid refinement with WENO approximations of the spatial derivative. Using 
the same problem as above we discretize the spatial grid with N = 101 equidistant points and use the ninth order 
weighted essentially non-oscillatory method (WEN09) [1] to differentiate the spatial derivatives. It is interesting 
to note that although the PDE we solve is linear, the use of the nonlinear method WEN09 results in a non-linear 


/K) = 
/'K) = 
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SSPRK 3,3 

2s3p 

2s4p 

3s5p 

A 

error 

Order 

error 

Order 

error 

Order 

error 

Order 

0.8 

7.99 X 10'® 

— 

1.86 X 10^® 

— 

1.96 X 10“® 

— 

6.47 X 10“* 

— 

0.7 

5.24 X 10“® 

3.16 

1.21 X 10“® 

3.17 

1.12 X 10“® 

4.17 

3.24 X 10“® 

5.17 

0.6 

3.27 X 10“® 

3.04 

7.61 X 10“® 

3.05 

6.02 X 10“’’ 

4.04 

1.49 X 10“® 

5.04 

0.5 

1.93 X lO”'^ 

2.87 

4.50 X 10“® 

2.88 

2.97 X 10“’’ 

3.87 

6.12 X 10“® 

4.87 

0.4 

9.70 X 10"® 

3.10 

2.25 X 10“® 

3.10 

1.18 X 10“’’ 

4.10 

1.96 X 10“® 

5.09 

0.3 

4.09 X 10"® 

3.00 

9.50 X 10“’’ 

3.00 

3.76 X 10“® 

3.99 

4.66 X 10““ 

5.00 

0.2 

1.21 X 10“® 

2.99 

2.81 X 10“’’ 

3.00 

7.43 X 10“® 

4.01 

6.13 X 10““ 

5.00 

0.1 

1.50 X 10"'^ 

3.01 

3.49 X 10“® 

3.01 

4.61 X 10““ 

4.00 

1.90 X 10““ 

5.01 

0.05 

1.88 X 10“® 

2.99 

4.36 X 10“® 

3.00 

2.88 X 10““ 

3.97 

5.97 X 10““ 

4.99 


Table 4: Convergence study for Example 3a, the linear advection problem with pseudospectral differentia¬ 
tion of the spatial derivatives. Here we use iV = 41 equidistant points between (0, 2tt), and At = XAx. The 
solution is evolved forward to time Tf = 2.0 using the explicit SSP Runge-Kutta method (SSPRK3,3), 
the two-stage third order two derivative method (2s3p), the two-stage fourth order two derivative method 
(2s4p), and three-stage fifth order two derivative method (3s5p). 


ODE. To evolve this ODE in time we use a range of time steps defined by At = AAa: for A = 0.9, 0.8, 0.7,0.6,0.5,0.4 
to compute the solution to final time Tf = 2.0 In Table 5 we list the errors for the SSPRK3,3, the two derivative 
two stage third order method in Section 3.2 with K ~ (listed as 2s3p in the table), the unique the two derivative 
two stage fourth order method (2s4p), and the the two derivative three stage fifth order method (22) (3s5p). We 
observe that the design-order of each method is verified and that once again the SSPRK3,3 method has larger errors 
than the multiderivative method 2s3p. These results indicate that the approximation of the second derivative term 
Ft via a Lax-Wendroff procedure and discretization in space does not affect the observed order of the time-stepping 
method, as long as the spatial errors do not dominate. 

To see what happens if the spatial discretization errors dominate over the time errors, we compare two different 
WENO spatial discretizations: the fifth order method WEN05 and the ninth order method WEN09. Here, we 
use N — 301 points in space, and we choose At = AAa; for A = 0.8, 0.6, 0.4,0.2,0.1,0.05. We evolve the solution in 
time to Tf = 2.0 using the 2s3p and SSPRK3,3 methods. The log-log graph of the errors vs. At is given in Figure 
8 We observe that the errors from both time discretizations with the WEN05 spatial discretization (dotted lines) 
have the correct orders when At is larger and the time errors dominate, but as At gets smaller the spatial errors 
dominate convergence is lost. On the other hand, for the range of At studied, the time-errors dominate over the 
spatial errors when using the highly accurate WEN09 (solid line) and so convergence is not lost. We note that 
to see this behavior, the spatial approximation must be sufficiently accurate. In this case this is attained by using 
301 points in space and a high order spatial discretization. We also studied this problem with a co-refinement of 
the spatial and temporal grids, where we use At = 0.8Aa: for Ax = and N — 41, 81,161,321. Figure 9 shows 
that while for larger At the errors for WEN05 are larger than for WEN09, and this is more pronounced for the 
multi-derivative methods than for the explicit Runge-Kutta, once the grid is sufficiently refined the temporal error 
dominates and we see the third order convergence in time. 

It is interesting to note that in all these studies the explicit SSP Runge-Kutta method SSPRK3,3 and the 
two-derivative method 2s3p behave similarly, indicating that the approximation of the second derivative does not 
play a role in the loss of accuracy. 

Example 4: co-refinement of the spatial and temporal grids on linear advection with WEN07 
Example 4a: In this example, we compare the errors and order of convergence of the multiderivative methods on 
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Figure 8: WEN05 vs WEN09 with grid re¬ 
finement in time but not space Example 3b. On 
the x-axis are logio{At) and on the y-axis are 
logio{error). The errors from both time dis¬ 
cretizations with the WEN05 spatial discretiza¬ 
tion (dotted lines) have the correct orders when 
At is larger and the time errors dominate, but 
as At gets smaller the spatial errors dominate 
convergence is lost. 



Figure 9: Example 3b: WEN05 vs WEN09 
with refinement in both space and time, At = 
0.8Ax . On the x-axis are logio{At) and on the 
y-axis are logio{error). Here we see that for 
larger At the errors for WEN05 are larger than 
for WEN09 (especially for the multiderivative 
method) but once the grid is sufficiently refined 
the temporal order dominates and third order 
convergence in time is observed. 


a linear and nonlinear problem. For the linear problem, we use the linear advection problem 

Ut + Ux = 0 on X G [—1,1] 

with the initial condition 

uo(x) = 0.5-I-0.5sin(7rx), x € [—1,1], (30) 

We compute both F and F using a seventh order WENO (WEN07) [1] spatial derivative. We set At = 0.8Ax 
and evolve the solution to final time of Tfi„ai = 2.0, at which point the exact solution is identical to the initial 
state. For time-stepping, we use the same three multiderivative schemes 2s3p, 2s4p, and 3s5p and the explicit SSP 
Runge-Kutta method SSPRK3,3. In Table 6, where we compare the errors and orders of these three methods. 
These numerical experiments show that once the mesh is sufficiently refined, the design-order of accuracy is reached. 
The results verify that the proposed schemes are genuinely high-order accurate despite the fact that the second 
derivative is not an exact second derivative of the method of lines formulation. 

Example 4b: Next, we present results for the more difficult non-linear Burgers equation, with initial conditions 
prescribed by 

17o(x) = 1.0-I-0.2sin(7rx), x € [—1,1], (31) 

and periodic boundary conditions. Once again, F and F are computed via the WEN07 spatial discretization. We 
use a constant time-step At — and run this problem with a final time of Tfi„ai = 1-4, before a shock 

forms. Since the solution at this point the solution remains smooth, we can use the method of characteristics to 
compute the the exact solution by 

U{t,x) = Uo{x-t-Uo{0), ^ = x-t-Uo{0. (32) 


20 











SSPRK 3,3 

2s3p 

2s4p 

3s5p 

A 

error 

Order 

error 

Order 

error 

Order 

error 

Order 

0.9 

7.37 X 10~® 

— 

1.71 X 10'® 

— 

8.25 X 10'® 

— 

1.24 X 10'® 

— 

0.8 

5.24 X 10“® 

2.89 

1.22 X 10'® 

2.89 

5.21 X 10'® 

3.89 

6.99 X 10'^® 

4.89 

0.7 

3.44 X 10“® 

3.14 

7.99 X 10'^ 

3.14 

3.00 X 10'® 

4.14 

3.52 X 10'^® 

5.14 

0.6 

2.18 X 10"® 

2.96 

5.08 X 10'^ 

2.96 

1.63 X 10'® 

3.96 

1.64 X 10'^® 

4.96 

0.5 

1.27 X 10"® 

2.98 

2.94 X 10'^ 

2.98 

7.89 X 10'® 

3.98 

6.61 X 10'“ 

4.97 

0.4 

6.47 X 10"'^ 

3.01 

1.50 X 10'^ 

3.01 

3.22 X 10'® 

4.01 

2.18 X 10'“ 

4.97 


Table 5: Convergence study for Example 3b, the linear advection problem with WEN09 differentiation of 
the spatial derivatives. Here we use N = 101 equidistant points between (0, 27 r), and At = AAx. 



SSPRK 3,3 

2s3p 

2s4p 

3s5p 

N 

error 

order 

error 

order 

error 

order 

error 

order 

41 

3.00 X 10'®'‘ 

— 

5.94 X 10'°® 

— 

7.54 X 10'°® 

— 

2.59 X 10'°® 

— 

81 

3.75 X 10'°® 

3.00 

7.39 X 10'°® 

3.01 

4.71 X 10'°’’ 

4.00 

8.03 X 10'°® 

5.01 

161 

4.69 X 10'°® 

3.00 

9.23 X 10'°^ 

3.00 

2.94 X 10'°® 

4.00 

2.51 X 10'°® 

5.00 

321 

5.86 X 10'°'^ 

3.00 

1.15 X 10'°^ 

3.00 

1.84 X 10'°® 

4.00 

7.82 X 10'“ 

5.00 

641 

7.32 X 10'°® 

3.00 

1.44 X 10'°® 

3.00 

1.15 X 10'’° 

4.00 

2.45 X 10'“ 

5.00 

1281 

9.15 X 10'°® 

3.00 

1.80 X 10'°® 

3.00 

7.19 X 10'“ 

4.00 

7.75 X 10'“ 

4.98 


Table 6: Convergence study for Example 4a: linear advection with WEN07 spatial differentiation and 
refinement of both spatial and temporal grids. Here we use N points in space and At = O.SAcc. The 
four methods used to evolve the solution to time Tf = 2.0 are the explicit SSP Runge-Kutta method 
SSPRK3,3 and the three multiderivative methods 2s3p, 2s4p, and 3s5p. Despite the fact that the second 
derivative F is not the exact second derivative Ft, we still observe high-order accuracy for all methods 
once the grids are sufficiently refined. 



SSPRK 3,3 

2s3p 

2s4p 

3s5p 

N 

error 

order 

error 

order 

error 

order 

error 

order 

161 

3.09 X 10'°^ 

— 

1.91 X lO'®"* 

— 

1.58 X 10'°'’ 

— 

1.67 X lO'®’ 

— 

321 

5.11 X 10'°® 

2.60 

2.00 X 10'°® 

3.25 

1.33 X 10'°® 

3.57 

1.76 X 10'°® 

3.25 

641 

6.96 X 10'°® 

2.88 

1.70 X 10'°® 

3.56 

7.12 X 10'°’’ 

4.22 

9.34 X lO'®’ 

4.23 

1281 

8.80 X 10'°’’ 

2.98 

1.81 X lO'®’ 

3.23 

3.38 X 10'°® 

4.39 

2.73 X 10'°® 

5.09 

2561 

1.10 X 10'°’’ 

3.00 

2.22 X 10'°® 

3.03 

1.99 X 10'°® 

4.09 

7.49 X 10'“ 

5.19 

5121 

1.38 X 10'°® 

3.00 

2.76 X 10'°® 

3.01 

1.24 X 10'“ 

4.01 

2.21 X 10'“ 

5.08 

10241 

1.72 X 10'°® 

3.00 

3.44 X 10'“ 

3.00 

7.75 X 10'“ 

4.00 

7.04 X 10'“ 

4.97 


Table 7: Convergence study for Example 4b: Burgers’ equation with WEN07 spatial differentiation and 
refinement of both spatial and temporal grids. Here we use N points in space and At = O.SAx. The 
four methods used to evolve the solution to time Tf = 2.0 are the explicit SSP Runge-Kutta method 
SSPRK3,3 and the three multiderivative methods 2s3p, 2s4p, and 3s5p. A more refined grid is needed in 
this example compared to the linear example, but for a sufficiently refined grid we still observe high-order 
accuracy for all methods. 
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We solve for the implicit variable ^ using Netwon iteration with a tolerance of 10“ The errors and order are 
presented in Table 7. We see that it takes an even smaller mesh size than the linear problem before the errors 
reach the asymptotic regime, but that once this happens we achieve the expected order. Again, these numerical 
experiments further validate the fact that if the spatial error is not allowed to dominate, the high-order design- 
accuracy of the time discretization attained despite the fact that we do not directly differentiate the method of 
lines formulation to define the second derivative. 


5 Conclusions 

With the increasing popularity of multi-stage multiderivative methods for use as time-stepping methods for hyper¬ 
bolic problems, the question of their strong stability properties needs to be addressed. In this work we presented 
an SSP formulation for multistage two-derivative methods. We assumed that, in addition to the forward Euler 
condition, the spatial discretization of interest satisfies a second derivative condition of the form (12). With these 
assumptions in mind, we formulated an optimization problem which enabled us to find optimal explicit SSP multi¬ 
stage two-derivative methods of up to order five, thus breaking the SSP order barrier for explicit SSP Runge-Kutta 
methods. Numerical test cases verify the convergence of these methods at the design-order, show that sharpness 
of the SSP condition in many cases, and demonstrate the need for SSP time-stepping methods in simulations 
where the spatial discretization is specially designed to satisfy certain nonlinear stability properties. Future work 
will involve building SSP multiderivative methods while assuming different base conditions (as in Remark 1) and 
with higher derivatives. Additional work will involve developing new spatial discretizations suited for use with 
SSP multiderivative time stepping methods. These methods will be based on WENO or discontinuous Galerkin 
methods and will satisfy pseudo-TVD and similar properties for systems of equations. 
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A Coefficients of three stage fourth order methods 


1. For A = I we obtain an SSP coefficient C = r — 1.1464. The Butcher array coefficients are given by 


0 

0 

0 


■ 0.528992280543542 ' 

0.436148675945340 

0 

0 

, 6 = 

0.105732787708912 

0.546571371212865 

0.156647174804152 

0 


_ 0.365274931747546 _ 

0 

0 

0 ■ 


0.074866026156687 ‘ 

0.095112833764436 

0.071032477596813 

0 

0.107904226252921 

-1 

o o 

, b = 

0.073410341982927 
0.048740310097159 _ 


The Shu-Osher arrays are Rb — (1,0,0,0)^, 


P = 


0 

0 

0 

0 

0.5000 

0 

0 

0 

0.253176729307242 

0.179580018745470 

0 

0 

0.181986129712082 

0.000000001889650 

0.418750476492704 

0 


Q 


0 

0 

0 

0 

0.5 

0 

0 

0 

0 

0.567243251947287 

0 

0 

0.140002406985210 

0.003037359947341 

0.256223624973014 

0 
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2. For K = 



we obtain an SSP coefficient C — r = 1.3927 Butcher formulation 


0 

0 

0 ■ 


' 0.515040964378407 ‘ 

0.443752012194422 

0 

0 

, b = 

0.178821699719783 

0.543193299768317 

0.149202742858795 

0 _ 


_ 0.306137335901811 _ 

0 

0 

0 1 


■ 0.072864982225864 ' 

0.098457924163299 

0.062758211639901 

0 

0.110738910914425 

-1 

o o 

, b = 

0.073840478463180 
_ 0.061973770357455 _ 


The Shu-Osher arrays are Rb — (1,0,0,0)^, 


0 

0 

0 

0.618033988749895 

0 

0 

0.362588515112176 

0.207801573327953 

0 

0.144580879241747 

0.110491604448675 

0.426371652664792 


0 ■ 
0 
0 
0 


■ 0 

0.381966011250105 

0 

0.078129569197367 


0 

0 

0.429609911559871 

0 


0 

0 

0 

0.240426294447419 


0 ■ 
0 
0 
0 


3. For if = 1 we obtain an SSP coefficient C = r — 1.6185. The Butcher array coefficients are given by 


0 

0 

0 ■ 


■ 0.502519798444212 ‘ 

0.452297224196082 

0 

0 

, b = 

0.210741084344740 

0.528050722182308 

0.159236998008155 

0 _ 


_ 0.286739117211047 _ 

0 

0 

0 1 


■ 0.071256397204544 ' 

0.102286389507741 

0.055482128781494 

0 

0.108677624192402 

-1 

o o 

, b = 

0.069475972085130 
_ 0.066877749079721 _ 


The Shu-Osher arrays are Rb — (1,0,0,0)^, 


0 

0 

0 

0 

0.732050807568877 

0 

0 

0 

0.457580533944888 

0.257727809835459 

0 

0 

0.137886171629970 

0.176326540063367 

0.464092174540814 

0 


0 

0 

0 

0.267949192431123 

0 

0 

0 

0.284691656219654 

0 

0.046502314960818 

0 

0.175192798805030 


0 ■ 
0 
0 
0 


B Two stage third order method 

This code gives the SSP coefficient and the Butcher and Shu Osher arrays for the optimal explicit SSP two stage 
third order method given the value K. 
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clear all 

k=sqrt(0.5); ’/, Choose K 
tab= [] ; 

y. Set up the polynomial for Q(3,l) to solve for SSP coefficient r 
AA=sqrt(k~2+2) - k; 
pO=2*k*(AA-2*k) + 4*k~3*AA; 
pl=-pO; 

p2=(l-p0)/(2*k~2); 

p3= -(p0/(2*k) + k)/(6*k~3); 

CC= [p3,p2 ,pl ,p0] ; "/, polynomial coefficients 
RC=roots(CC); 

r=RC(f indCabs (imag(RC)) <10~-15) ) ; "/, SSP coefficient is the only real root. 

% - 

y. Once we have the k and r we want we define the method following (21) 
a= (k*sqrt(k~2+2)-k~2)/r; 

b2 = ((k-2*(l-l/r)) + r*(.5-l/(6*a)))/(k-2+.5*r*a); 

bl=l-b2; 

ahat=.5*a~2; 

bhatl=.5*(l-b2*a)-1/(6*a) ; 
bhat2=l/(6*a)-.5*b2*a; 
y. The Butcher arrays are given by 
A=zeros(2,2); Ahat=A; 

A(2,l)=a; 
b=[bl,b2] ; 

Ahat(2,1)=ahat; 
bhat=[bhatl,bhat2]; 

S=[0 0 0 ; a 0 0; bl b2 0] ; 

Shat= [0 0 0 ; ahat 0 0 ; bhatl bhat2 0] ; 
y. The Shu-Osher matrices are given by 
I=eye(3);e=ones(3,1); 

Ri=I+r*S+(r~2/k~2)*Shat; 
v=Ri\e; 

P = r*(Ri\S) ; 
q= r-2/k-2*(Ri\Shat); 

violation=min(min( [v, S, Shat, P, Q])) "/, use this to check that all these are positive 
tab= [tab; [k,r, violation] ] ; "/ build the table of values 


C Three stage fifth order method 

This code gives the SSP coefficient and the Butcher and Shu Osher arrays for the optimal explicit SSP three stage 
fifth order method given the value K. 

clear all 
format long 
syms r 
tab= [] ; 

k=sqrt(0.5) "/The second derivative condition coefficient K 
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°/« Find the SSP coefficient C given K 

a21= 240*k-6*(l -r - r-2/(2*k~2) + r-3/(6*k-2) + r-4/(24*k-4) - r~5/(120*k-4))/r-6; 

Q31=10*r"2*a21~4 - 100*k~2*a21~3 - 10*r"2*a21~3 + 130*k~2*a21~2 + 3*r~2*a21~2 - 50*k"2*a21 +6*k~2; 

RC=vpasolve(simplify(r~22*Q31)==0); 
rr=RC(find(abs(imag(RC))<10~-15)); 

C= max(rr) "/,The SSP coefficient 

% - 

y. The Butcher array coefficients given K and C 

a21= 240*k-6*(l -C - C-2/(2*k~2) + C-3/(6*k-2) + C-4/(24*k-4) - C~5/(120*k-4))/C-6; 

ah32= ( (3/5 -a21)“2/(a21*(l-2*a21)“3) - (3/5 -a21)/(l-2*a21)"2 )/10; 

ah31= ( (3/5 -a21)-2/(l-2*a21)-2)/2 -ah32; 

a31= (3/5 -a21)/(l-2*a21); 

bh2=(2*a31-l)/(12*a21*(a31-a21)); 

bh3=(l-2*a21)/(12*a31*(a31-a21)); 

bhl=l/2-bh2-bh3; 

ah21=(l/24 - bh3*(ah31+ah32) )/bh2; 
y. Build the Butcher matrices 
a= C*a21+ah21*C-2/k~2; 
b= C*a31+ah31*C-2/k~2; 
c= ah32*C~2/k"2; 
d= C+ bhl*C~2/k-2; 
e= bh2*C-2/k-2; 
f = bh3*C~2/k-2; 

Ri=([l 00 0;al00;bcl0;def 1]); 

S=[0 0 0 0; a21 0 0 0; a31 0 0 0 ; 1 0 0 0] ; 

Shat=[0 000; ah21 000; ah31 ah32 0 0 ;bhl bh2 bh3 0]; 
y. The Shu Qsher matrices given C 
eone=ones(4,1) 
v=Ri\eone; 

P = C*(Ri\S); 
q= C-2/k-2*(Ri\Shat) ; 

y. Double check that there are no violations of the SSP conditions: 

violation=min(min( [v, S, Shat, P, Q])) "/, use this to check that all these are positive 
tab= [tab; [k,a21 ,C,violation]] ; "/, build the table of values 
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